#include "gmshImport.hpp"
#include "gmshImport_testMSH.hpp"
#include "util.hpp"

#include <catch2/catch_test_macros.hpp>
#include <catch2/catch_approx.hpp>
#include <sstream>

using namespace hms;

using Catch::Approx;

TEST_CASE( "readFile" ){
	std::stringstream ss;
	
	ss << testMSH;
	gmsh::MSH msh { gmsh::readFile( ss ) };
	
	SECTION("Mesh format"){
		REQUIRE( msh.version == Approx(4.1) );
		REQUIRE( !msh.isBinary );
		REQUIRE( msh.dataSize == 8 );
	}
	SECTION("Physical names"){
		REQUIRE( msh.physicalNames.size() == 4 );
		REQUIRE( msh.physicalNames[0].dim == gmsh::Dim::curve );
		REQUIRE( msh.physicalNames[1].dim == gmsh::Dim::curve );
		REQUIRE( msh.physicalNames[2].dim == gmsh::Dim::curve );
		REQUIRE( msh.physicalNames[3].dim == gmsh::Dim::curve );
		REQUIRE( msh.physicalNames[0].tag == 1 );
		REQUIRE( msh.physicalNames[1].tag == 2 );
		REQUIRE( msh.physicalNames[2].tag == 3 );
		REQUIRE( msh.physicalNames[3].tag == 4 );
		REQUIRE( msh.physicalNames[0].name == "bottom" );
		REQUIRE( msh.physicalNames[1].name == "right" );
		REQUIRE( msh.physicalNames[2].name == "top" );
		REQUIRE( msh.physicalNames[3].name == "left" );
	}
	SECTION("Entities"){
		SECTION("Points"){
			REQUIRE( msh.points.size() == 4 );
			REQUIRE( msh.points[0].tag == 1 );
			REQUIRE( msh.points[1].tag == 2 );
			REQUIRE( msh.points[2].tag == 3 );
			REQUIRE( msh.points[3].tag == 4 );
			REQUIRE( msh.points[0].physicalTags.empty() );
			REQUIRE( msh.points[1].physicalTags.empty() );
			REQUIRE( msh.points[2].physicalTags.empty() );
			REQUIRE( msh.points[3].physicalTags.empty() );
			REQUIRE( msh.points[0].coords == Vector3s{0  , 0  , 0} );
			REQUIRE( msh.points[1].coords == Vector3s{0.6, 0  , 0} );
			REQUIRE( msh.points[2].coords == Vector3s{0.6, 0.5, 0} );
			REQUIRE( msh.points[3].coords == Vector3s{0  , 0.5, 0} );
		}
		SECTION("Curves"){
			REQUIRE( msh.curves.size() == 4 );
			REQUIRE( msh.curves[0].tag == 1 );
			REQUIRE( msh.curves[1].tag == 2 );
			REQUIRE( msh.curves[2].tag == 3 );
			REQUIRE( msh.curves[3].tag == 4 );
			REQUIRE( msh.curves[0].physicalTags == IdxVector{1} );
			REQUIRE( msh.curves[1].physicalTags == IdxVector{2} );
			REQUIRE( msh.curves[2].physicalTags == IdxVector{3} );
			REQUIRE( msh.curves[3].physicalTags == IdxVector{4} );
			REQUIRE( msh.curves[0].minCoords == Vector3s{0  , 0  , 0} );
			REQUIRE( msh.curves[1].minCoords == Vector3s{0.6, 0  , 0} );
			REQUIRE( msh.curves[2].minCoords == Vector3s{0  , 0.5, 0} );
			REQUIRE( msh.curves[3].minCoords == Vector3s{0  , 0  , 0} );
			REQUIRE( msh.curves[0].maxCoords == Vector3s{0.6, 0  , 0} );
			REQUIRE( msh.curves[1].maxCoords == Vector3s{0.6, 0.5, 0} );
			REQUIRE( msh.curves[2].maxCoords == Vector3s{0.6, 0.5, 0} );
			REQUIRE( msh.curves[3].maxCoords == Vector3s{0  , 0.5, 0} );
			REQUIRE( msh.curves[0].pointTags == IdxVector{1, -2} );
			REQUIRE( msh.curves[1].pointTags == IdxVector{2, -3} );
			REQUIRE( msh.curves[2].pointTags == IdxVector{3, -4} );
			REQUIRE( msh.curves[3].pointTags == IdxVector{4, -1} );
		}
	}
	SECTION("Nodes"){
		REQUIRE( msh.nNodes == 30 );
		REQUIRE( msh.nodeBlocks.size() == 9 );
		REQUIRE( msh.nodeBlocks[0].entityDim == gmsh::Dim::point );
		REQUIRE( msh.nodeBlocks[1].entityDim == gmsh::Dim::point );
		REQUIRE( msh.nodeBlocks[2].entityDim == gmsh::Dim::point );
		REQUIRE( msh.nodeBlocks[3].entityDim == gmsh::Dim::point );
		REQUIRE( msh.nodeBlocks[4].entityDim == gmsh::Dim::curve );
		REQUIRE( msh.nodeBlocks[5].entityDim == gmsh::Dim::curve );
		REQUIRE( msh.nodeBlocks[6].entityDim == gmsh::Dim::curve );
		REQUIRE( msh.nodeBlocks[7].entityDim == gmsh::Dim::curve );
		REQUIRE( msh.nodeBlocks[8].entityDim == gmsh::Dim::surface );
		REQUIRE( msh.nodeBlocks[0].entityTag == 1 );
		REQUIRE( msh.nodeBlocks[1].entityTag == 2 );
		REQUIRE( msh.nodeBlocks[2].entityTag == 3 );
		REQUIRE( msh.nodeBlocks[3].entityTag == 4 );
		REQUIRE( msh.nodeBlocks[4].entityTag == 1 );
		REQUIRE( msh.nodeBlocks[5].entityTag == 2 );
		REQUIRE( msh.nodeBlocks[6].entityTag == 3 );
		REQUIRE( msh.nodeBlocks[7].entityTag == 4 );
		REQUIRE( msh.nodeBlocks[8].entityTag == 6 );
		REQUIRE( msh.nodeBlocks[0].nodeTags == IdxVector{1} );
		REQUIRE( msh.nodeBlocks[1].nodeTags == IdxVector{2} );
		REQUIRE( msh.nodeBlocks[2].nodeTags == IdxVector{3} );
		REQUIRE( msh.nodeBlocks[3].nodeTags == IdxVector{4} );
		REQUIRE( msh.nodeBlocks[4].nodeTags == IdxVector{ 5,  6,  7,  8} );
		REQUIRE( msh.nodeBlocks[5].nodeTags == IdxVector{ 9, 10, 11    } );
		REQUIRE( msh.nodeBlocks[6].nodeTags == IdxVector{12, 13, 14, 15} );
		REQUIRE( msh.nodeBlocks[7].nodeTags == IdxVector{16, 17, 18    } );
		REQUIRE( msh.nodeBlocks[8].nodeTags ==
			IdxVector{19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30} );
		Matrix3Xs coords (3, 30);
		coords <<
			0,0.6,0.6,0  ,0.12,0.24,0.36,0.48,0.6  ,0.6 ,0.6  ,0.48,0.36,0.24,0.12,0    ,0   ,0    ,0.12 ,0.12,0.12 ,0.24 ,0.24,0.24 ,0.36 ,0.36,0.36 ,0.48 ,0.48,0.48 ,
			0,0  ,0.5,0.5,0   ,0   ,0   ,0   ,0.125,0.25,0.375,0.5 ,0.5 ,0.5 ,0.5 ,0.375,0.25,0.125,0.125,0.25,0.375,0.125,0.25,0.375,0.125,0.25,0.375,0.125,0.25,0.375,
			0,0  ,0  ,0  ,0   ,0   ,0   ,0   ,0    ,0   ,0    ,0   ,0   ,0   ,0   ,0    ,0   ,0    ,0    ,0   ,0    ,0    ,0   ,0    ,0    ,0   ,0    ,0    ,0   ,0    ;
		REQUIRE( msh.nodeBlocks[0].coords.isApprox(coords.col(0), 1e-6 ) );
		REQUIRE( msh.nodeBlocks[1].coords.isApprox(coords.col(1), 1e-6 ) );
		REQUIRE( msh.nodeBlocks[2].coords.isApprox(coords.col(2), 1e-6 ) );
		REQUIRE( msh.nodeBlocks[3].coords.isApprox(coords.col(3), 1e-6 ) );
		REQUIRE( msh.nodeBlocks[4].coords.isApprox(coords.block(0, 4,3, 4), 1e-6) );
		REQUIRE( msh.nodeBlocks[5].coords.isApprox(coords.block(0, 8,3, 3), 1e-6) );
		REQUIRE( msh.nodeBlocks[6].coords.isApprox(coords.block(0,11,3, 4), 1e-6) );
		REQUIRE( msh.nodeBlocks[7].coords.isApprox(coords.block(0,15,3, 3), 1e-6) );
		REQUIRE( msh.nodeBlocks[8].coords.isApprox(coords.block(0,18,3,12), 1e-6) );
	}
	SECTION("Elements"){
		REQUIRE( msh.elementBlocks.size() == 9 );
		REQUIRE( msh.elementBlocks[0].entityDim == gmsh::Dim::point );
		REQUIRE( msh.elementBlocks[1].entityDim == gmsh::Dim::point );
		REQUIRE( msh.elementBlocks[2].entityDim == gmsh::Dim::point );
		REQUIRE( msh.elementBlocks[3].entityDim == gmsh::Dim::point );
		REQUIRE( msh.elementBlocks[4].entityDim == gmsh::Dim::curve );
		REQUIRE( msh.elementBlocks[5].entityDim == gmsh::Dim::curve );
		REQUIRE( msh.elementBlocks[6].entityDim == gmsh::Dim::curve );
		REQUIRE( msh.elementBlocks[7].entityDim == gmsh::Dim::curve );
		REQUIRE( msh.elementBlocks[8].entityDim == gmsh::Dim::surface );
		REQUIRE( msh.elementBlocks[0].entityTag == 1 );
		REQUIRE( msh.elementBlocks[1].entityTag == 2 );
		REQUIRE( msh.elementBlocks[2].entityTag == 3 );
		REQUIRE( msh.elementBlocks[3].entityTag == 4 );
		REQUIRE( msh.elementBlocks[4].entityTag == 1 );
		REQUIRE( msh.elementBlocks[5].entityTag == 2 );
		REQUIRE( msh.elementBlocks[6].entityTag == 3 );
		REQUIRE( msh.elementBlocks[7].entityTag == 4 );
		REQUIRE( msh.elementBlocks[8].entityTag == 6 );
		REQUIRE( msh.elementBlocks[0].elementType == gmsh::ElementType::point );
		REQUIRE( msh.elementBlocks[1].elementType == gmsh::ElementType::point );
		REQUIRE( msh.elementBlocks[2].elementType == gmsh::ElementType::point );
		REQUIRE( msh.elementBlocks[3].elementType == gmsh::ElementType::point );
		REQUIRE( msh.elementBlocks[4].elementType == gmsh::ElementType::line  );
		REQUIRE( msh.elementBlocks[5].elementType == gmsh::ElementType::line  );
		REQUIRE( msh.elementBlocks[6].elementType == gmsh::ElementType::line  );
		REQUIRE( msh.elementBlocks[7].elementType == gmsh::ElementType::line  );
		REQUIRE( msh.elementBlocks[8].elementType == gmsh::ElementType::quadrilateral);
		REQUIRE( msh.elementBlocks[0].elementTags == IdxVector{1} );
		REQUIRE( msh.elementBlocks[1].elementTags == IdxVector{2} );
		REQUIRE( msh.elementBlocks[2].elementTags == IdxVector{3} );
		REQUIRE( msh.elementBlocks[3].elementTags == IdxVector{4} );
		REQUIRE( msh.elementBlocks[4].elementTags == IdxVector{ 5, 6, 7, 8, 9});
		REQUIRE( msh.elementBlocks[5].elementTags == IdxVector{10,11,12,13   });
		REQUIRE( msh.elementBlocks[6].elementTags == IdxVector{14,15,16,17,18});
		REQUIRE( msh.elementBlocks[7].elementTags == IdxVector{19,20,21,22   });
		REQUIRE( msh.elementBlocks[8].elementTags == IdxVector{
			23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42
		});
		ArrayXXi nodeTags (1,4);
		nodeTags << 1,2,3,4;
		REQUIRE( ( msh.elementBlocks[0].nodeTags == nodeTags.col(0) ).all() );
		REQUIRE( ( msh.elementBlocks[1].nodeTags == nodeTags.col(1) ).all() );
		REQUIRE( ( msh.elementBlocks[2].nodeTags == nodeTags.col(2) ).all() );
		REQUIRE( ( msh.elementBlocks[3].nodeTags == nodeTags.col(3) ).all() );
		nodeTags.resize(2,5);
		nodeTags <<
			1, 5, 6, 7, 8,
			5, 6, 7, 8, 2;
		REQUIRE( ( msh.elementBlocks[4].nodeTags == nodeTags ).all() );
		nodeTags.resize(2,4);
		nodeTags <<
			2,  9, 10, 11,
			9, 10, 11,  3;
		REQUIRE( ( msh.elementBlocks[5].nodeTags == nodeTags ).all() );
		nodeTags.resize(2,5);
		nodeTags <<
			 3, 12, 13, 14, 15,
			12, 13, 14, 15,  4; 
		REQUIRE( ( msh.elementBlocks[6].nodeTags == nodeTags ).all() );
		nodeTags.resize(2,4);
		nodeTags <<
			 4, 16, 17, 18,
			16, 17, 18,  1;
		REQUIRE( ( msh.elementBlocks[7].nodeTags == nodeTags ).all() );
		nodeTags.resize(4,20);
		nodeTags <<
			 1,18,17,16, 5,19,20,21, 6,22,23,24, 7,25,26,27, 8,28,29,30,
			 5,19,20,21, 6,22,23,24, 7,25,26,27, 8,28,29,30, 2, 9,10,11,
			19,20,21,15,22,23,24,14,25,26,27,13,28,29,30,12, 9,10,11, 3,
			18,17,16, 4,19,20,21,15,22,23,24,14,25,26,27,13,28,29,30,12;
		REQUIRE( ( msh.elementBlocks[8].nodeTags == nodeTags ).all() );
	}
}

TEST_CASE( "gmshImport, quadrilaterals" ){
	using hms::exactlyEqual;
	
	std::stringstream ss;
	ss << testMSH;
	
	hms::Mesh mesh { gmsh::gmshToHms( ss ) };
	
	SECTION( "Info" ){
		REQUIRE( mesh.nDomainCells() == 20 );
		REQUIRE( mesh.nCells() == 38 );
		REQUIRE( mesh.nVertices() == 30 );
		REQUIRE( mesh.nEdges() == 49 );
		REQUIRE( mesh.nInnerEdges() == 31 );
		REQUIRE( mesh.nGhostCells() == 18 );
		REQUIRE( mesh.domainArea() == Approx(0.3) );
	}

	// Array4Xi domainCellVerticesStored ( 4, mesh.nDomainCells() );
	// ArrayNXs<8> domainCellVertexCoords ( 8, mesh.nDomainCells() );
	// for ( int i=0; i<mesh.nDomainCells(); ++i ){
	// 	domainCellVerticesStored.col(i) = mesh.cell(i).vertices(true);
	// 	domainCellVertexCoords.col(i) = mesh.vertices()(
	// 		hms::all, domainCellVerticesStored.col(i)
	// 	).reshaped();
	// }

	// std::cout << "Vertices:\n" << mesh.vertices()
	// 	<< "\nEdge vertex indices:\n" << mesh.edges().vertices()
	// 	<< "\nEdge lengths:\n" << mesh.edges().lengths().transpose()
	// 	<< "\nEdge adjacent cells:\n" << mesh.edges().cells()
	// 	<< "\nCell vertex indices:\n" << domainCellVerticesStored
	// 	<< "\nCell vertex coords:\n" << domainCellVertexCoords
	// 	<< "\nCell edge indices:\n" << mesh.cells().domainCellEdges()
	// 	<< "\nCell areas:\n" << mesh.cells().areas().transpose()
	// 	<< "\nCell centroids:\n" << mesh.cells().centroids()
	// 	<< "\nCell circumfs:\n" << mesh.cells().circumfs().transpose()
	// 	<< "\n";

	auto check = [](
		const auto& target,
		const auto& actual
	){
		bool allFound {true};
		for ( const auto& edge : actual.colwise() ){
			bool found {false};
			for ( const auto& candidate : target.colwise() ){
				if ( hms::sameEntries(edge, candidate) ){
					found = true;
					break;
				}
			}
			if ( !found ){
				allFound = false;
				break;
			}
		}
		REQUIRE( allFound );
	};

	SECTION( "Edge vertices" ){

		Array2Xi vertsTarget { ArrayNi<98>{
			5, 19,
			19, 18,
			19, 20,
			20, 17,
			20, 21,
			21, 16,
			21, 15,
			6, 22,
			22, 19,
			22, 23,
			23, 20,
			23, 24,
			24, 21,
			24, 14,
			7, 25,
			25, 22,
			25, 26,
			26, 23,
			26, 27,
			27, 24,
			27, 13,
			8, 28,
			28, 25,
			28, 29,
			29, 26,
			29, 30,
			30, 27,
			30, 12,
			9, 28,
			10, 29,
			11, 30,

			1,  5,
			5,  6,
			6,  7,
			7,  8,
			8,  2,
			2,  9,
			9, 10,
			10, 11,
			11,  3,
			3, 12,
			12, 13,
			13, 14,
			14, 15,
			15,  4,
			4, 16,
			16, 17,
			17, 18,
			18,  1,
		}.reshaped(2, 49) };
		/* the +1 is because tags in MSH files are 1-based, whereas here,
		 * indices are 0-based */
		vertsTarget -= 1;

		SECTION( "Inner edges" ){
			check(
				vertsTarget.leftCols( mesh.nInnerEdges() ),
				mesh.edges().vertices().leftCols( mesh.nInnerEdges() )
			);
		}
		SECTION( "Boundary edges" ){
			check(
				vertsTarget.rightCols( mesh.nBoundaryEdges() ),
				mesh.edges().vertices().rightCols( mesh.nBoundaryEdges() )
			);
		}
	}
	SECTION( "Edge parent cells" ){
		Array2Xi cellsTarget { ArrayNi<98>{
			23,27,
			23,24,
			24,28,
			24,25,
			25,29,
			25,26,
			26,30,
			27,31,
			27,28,
			28,32,
			28,29,
			29,33,
			29,30,
			30,34,
			31,35,
			31,32,
			32,36,
			32,33,
			33,37,
			33,34,
			34,38,
			35,39,
			35,36,
			36,40,
			36,37,
			37,41,
			37,38,
			38,42,
			39,40,
			40,41,
			41,42,

			23,43,
			27,44,
			31,45,
			35,46,
			39,47,
			39,48,
			40,49,
			41,50,
			42,51,
			42,52,
			38,53,
			34,54,
			30,55,
			26,56,
			26,57,
			25,58,
			24,59,
			23,60,
		}.reshaped(2,49) };
		/* in gmsh, elements share the same pool of tags, so that cell indices
		 * start with the next tag after the highest boundary line tag. With
		 * four points and 18 lines (->22), the next tag is 23 */
		cellsTarget -= 23;

		SECTION( "Inner edges" ){
			check(
				cellsTarget.leftCols( mesh.nInnerEdges() ),
				mesh.edges().cells().leftCols( mesh.nInnerEdges() )
			);
		}
		SECTION( "Boundary edges" ){
			// check(
			// 	cellsTarget.rightCols( mesh.nBoundaryEdges() ),
			// 	mesh.edges().cells().rightCols( mesh.nBoundaryEdges() )
			// );
			for (int i=0; i<2; ++i)
				REQUIRE( sameEntries(
					mesh.edges().cells().rightCols( mesh.nBoundaryEdges() ).row(i),
					cellsTarget.rightCols( mesh.nBoundaryEdges() ).row(i)
				) );
		}
	}
	SECTION( "cell vertices" ){
		Array4Xi domainVerticesTarget { ArrayNi<4*20>{
			1, 5,19,18,
			18,19,20,17,
			17,20,21,16,
			16,21,15, 4,
			5, 6,22,19,
			19,22,23,20,
			20,23,24,21,
			21,24,14,15,
			6, 7,25,22,
			22,25,26,23,
			23,26,27,24,
			24,27,13,14,
			7, 8,28,25,
			25,28,29,26,
			26,29,30,27,
			27,30,12,13,
			8, 2, 9,28,
			28, 9,10,29,
			29,10,11,30,
			30,11, 3,12,
		}.reshaped(4,20) };
		Array2Xi ghostVerticesTarget { ArrayNi<18*2>{
			1, 5,
			5, 6,
			6, 7,
			7, 8,
			8, 2,
			2, 9,
			9,10,
			10,11,
			11, 3,
			12, 3,
			13,12,
			14,13,
			15,14,
			4,15,
			16, 4,
			17,16,
			18,17,
			1,18,
		}.reshaped(2,18) };
		domainVerticesTarget -= 1;
		ghostVerticesTarget -= 1;

		Array4Xi domainVertices (4, mesh.nDomainCells() );
		for ( int i=0; i<mesh.nDomainCells(); ++i )
			domainVertices.col(i) = mesh.cell(i).vertices(true);
		Array2Xi ghostVertices (2, mesh.nGhostCells() );
		for (int i=0; i<mesh.nGhostCells(); ++i)
			ghostVertices.col(i) =
				mesh.cell( i + mesh.nDomainCells() ).vertices(true);

		SECTION( "domain cells" ){
			check( domainVerticesTarget, domainVertices );
		}
		SECTION( "ghost cells" ){
			check( ghostVerticesTarget, ghostVertices );
		}
	}
	SECTION( "boundary patches" ){
		ArrayNi<5> bottom {  5,  6,  7,  8,  9 };
		ArrayNi<4> right  { 10, 11, 12, 13     };
		ArrayNi<5> top    { 14, 15, 16, 17, 18 };
		ArrayNi<4> left   { 19, 20, 21, 22     };
		bottom.array() += -5+31;
		right .array() += -5+31;
		top   .array() += -5+31;
		left  .array() += -5+31;
		REQUIRE( exactlyEqual(mesh.boundaryPatch("bottom").edges(), bottom ) );
		REQUIRE( exactlyEqual(mesh.boundaryPatch("right" ).edges(), right  ) );
		REQUIRE( exactlyEqual(mesh.boundaryPatch("top"   ).edges(), top    ) );
		REQUIRE( exactlyEqual(mesh.boundaryPatch("left"  ).edges(), left   ) );
	}
}

TEST_CASE( "gmshImport, mixed" ){

	using hms::exactlyEqual;
	
	std::stringstream ss;
	ss << testMSH2;
	hms::Mesh mesh { gmsh::gmshToHms(ss) };
	
	SECTION( "Info" ){
		REQUIRE( mesh.nDomainCells() == 25 );
		REQUIRE( mesh.nCells() == 25+18 );
		REQUIRE( mesh.nVertices() == 18+17 );
		REQUIRE( mesh.nEdges() == 41+18 );
		REQUIRE( mesh.nInnerEdges() == 41 );
		REQUIRE( mesh.nGhostCells() == 18 );
		REQUIRE( mesh.domainArea() == Approx(0.3) );
	}

	auto check = [](
		const auto& target,
		const auto& actual
	){
		bool allFound {true};
		for ( const auto& edge : actual.colwise() ){
			bool found {false};
			for ( const auto& candidate : target.colwise() ){
				if ( hms::sameEntries(edge, candidate) ){
					found = true;
					break;
				}
			}
			if ( !found ){
				allFound = false;
				break;
			}
		}
		REQUIRE( allFound );
	};

	SECTION( "Edge vertices" ){

		Array2Xi vertsTarget { ArrayNi<2*59>{
			22,23,
			23,19,
			19,28,
			28,22,
			21,24,
			24,19,
			19,27,
			27,21,
			20,27,
			23,20,
			33,28,
			24,33,
			30,16,
			17,25,
			25,30,
			29, 9,
			10,26,
			26,29,
			32,25,
			18,32,
			31,26,
			11,31,
			23, 7,
			7,34,
			34,20,
			22, 6,
			14,33,
			24,13,
			21,35,
			35,13,
			12,31,
			30,15,
			29, 8,
			5,32,
			32,22,
			28,25,
			20,26,
			26,21,
			33,30,
			29,34,
			35,31,
			1, 5,
			5, 6,
			6, 7,
			7, 8,
			8, 2,
			2, 9,
			9,10,
			10,11,
			11, 3,
			3,12,
			12,13,
			13,14,
			14,15,
			15, 4,
			4,16,
			16,17,
			17,18,
			18, 1,
		}.reshaped(2, 59) };
		/* the +1 is because tags in MSH files are 1-based, whereas here,
		 * indices are 0-based */
		vertsTarget -= 1;

		SECTION( "Inner edges" ){
			check(
				vertsTarget.leftCols( mesh.nInnerEdges() ),
				mesh.edges().vertices().leftCols( mesh.nInnerEdges() )
			);
		}
		SECTION( "Boundary edges" ){
			check(
				vertsTarget.rightCols( mesh.nBoundaryEdges() ),
				mesh.edges().vertices().rightCols( mesh.nBoundaryEdges() )
			);
		}
	}
	SECTION( "Edge parent cells" ){
	/* 4 corners and 18 boundary lines take up the first 22 element tags, that's
	 * why 23 is added to the stored indices, to be better able to compare it
	 * when viewing the file in gmsh */
		Array2Xi cellsTarget { ArrayNi<2*59>{
			23,32,
			23,25,
			23,26,
			23,39,
			24,34,
			24,26,
			24,25,
			24,40,
			25,40,
			25,31,
			26,41,
			26,33,
			27,36,
			27,29,
			27,41,
			28,37,
			28,30,
			28,46,
			29,39,
			29,38,
			30,47,
			30,35,
			31,32,
			31,43,
			31,46,
			32,44,
			33,42,
			33,34,
			34,47,
			34,45,
			35,45,
			36,42,
			37,43,
			38,44,
			39,44,
			39,41,
			40,46,
			40,47,
			41,42,
			43,46,
			45,47,

			38,48,
			44,49,
			32,50,
			43,51,
			37,52,
			37,53,
			28,54,
			30,55,
			35,56,
			35,57,
			45,58,
			33,59,
			42,60,
			36,61,
			36,62,
			27,63,
			29,64,
			38,65,
		}.reshaped(2,59) };
		cellsTarget -= 23;

		SECTION( "Inner edges" ){
			check(
				cellsTarget.leftCols( mesh.nInnerEdges() ),
				mesh.edges().cells().leftCols( mesh.nInnerEdges() )
			);
		}
		SECTION( "Boundary edges" ){
			// check(
			// 	cellsTarget.rightCols( mesh.nBoundaryEdges() ),
			// 	mesh.edges().cells().rightCols( mesh.nBoundaryEdges() )
			// );
			for (int i=0; i<2; ++i)
				REQUIRE( sameEntries(
					mesh.edges().cells().rightCols( mesh.nBoundaryEdges() ).row(i),
					cellsTarget.rightCols( mesh.nBoundaryEdges() ).row(i)
				) );
		}
	}
	SECTION( "cell vertices" ){
		Array4Xi domainVerticesTarget { ArrayNi<4*25>{
			22,23,19,28,
			19,27,21,24,
			23,20,27,19,
			28,19,24,33,
			17,25,30,16,
			29, 9,10,26,
			32,25,17,18,
			10,11,31,26,
			7,34,20,23,
			6, 7,23,22,
			24,13,14,33,
			21,35,13,24,
			11, 3,12,31,
			16,30,15, 4,
			8, 2, 9,29,
			1, 5,32,18,
			32,22,28,25,
			20,26,21,27,
			25,28,33,30,
			33,14,15,30,
			7, 8,29,34,
			5, 6,22,32,
			35,31,12,13,
			34,29,26,20,
			26,31,35,21,
		}.reshaped(4,25) };
		Array2Xi ghostVerticesTarget { ArrayNi<18*2>{
			1, 5,
			5, 6,
			6, 7,
			7, 8,
			8, 2,
			2, 9,
			9,10,
			10,11,
			11, 3,
			12, 3,
			13,12,
			14,13,
			15,14,
			4,15,
			16, 4,
			17,16,
			18,17,
			1,18,
		}.reshaped(2,18) };
		/* account for 1-based vertices */
		domainVerticesTarget -= 1;
		ghostVerticesTarget -= 1;

		Array4Xi domainVertices (4, mesh.nDomainCells() );
		for ( int i=0; i<mesh.nDomainCells(); ++i )
			domainVertices.col(i) = mesh.cell(i).vertices(true);
		Array2Xi ghostVertices (2, mesh.nGhostCells() );
		for (int i=0; i<mesh.nGhostCells(); ++i)
			ghostVertices.col(i) =
				mesh.cell( i + mesh.nDomainCells() ).vertices(true);

		SECTION( "domain cells" ){
			check( domainVerticesTarget, domainVertices );
		}
		SECTION( "ghost cells" ){
			check( ghostVerticesTarget, ghostVertices );
		}
	}
	SECTION( "boundary patches" ){
		ArrayNi<5> bottom {  5,  6,  7,  8,  9 };
		ArrayNi<4> right  { 10, 11, 12, 13     };
		ArrayNi<5> top    { 14, 15, 16, 17, 18 };
		ArrayNi<4> left   { 19, 20, 21, 22     };
		bottom.array() += -5+41;
		right .array() += -5+41;
		top   .array() += -5+41;
		left  .array() += -5+41;
		REQUIRE( exactlyEqual(mesh.boundaryPatch("bottom").edges(), bottom ) );
		REQUIRE( exactlyEqual(mesh.boundaryPatch("right" ).edges(), right  ) );
		REQUIRE( exactlyEqual(mesh.boundaryPatch("top"   ).edges(), top    ) );
		REQUIRE( exactlyEqual(mesh.boundaryPatch("left"  ).edges(), left   ) );
	}
}

// /* usage */
// #include "meshTypes.hpp"
// #include "vtkMeshExport.hpp"
// #include <iostream>
// #include <filesystem>
// namespace fs = std::filesystem;

// TEST_CASE( "Importing gmsh file" ){

// fs::path dir { "/home/lenni/tubCloud/coding/gmsh-test" };
// std::ifstream file { dir / "trivialMeshTri.msh" };
// std::cout << "File is open: " << std::boolalpha << file.is_open() << "\n";
// hms::Mesh mesh { gmsh::gmshToHms(file) };
// mesh.printInfo();
// }
// TEST_CASE( "Importing gmsh file 2" ){
// 	using namespace hms;
// 	MeshVariant meshVar;

// 	fs::path dir { "/home/lenni/tubCloud/coding/hms_0_01/tests" };
// 	std::ifstream file { dir / "gmshImport_testMSH2.msh" };
// 	std::cout << "File is open: " << std::boolalpha << file.is_open() << "\n";
	
// 	hms::Mesh mesh { gmsh::gmshToHms(file) };
// 	mesh.printInfo();

// 	meshVar = std::move(mesh);

// 	vtk::meshToVtk(std::cout, meshVar);
// }